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A REVIEW  OF  MODERN  DEVELOPMENTS  IN  LARGE  EDDY 
SIMULATION  OF  TURBULENT  REACTIVE  FLOWS 


PEYMAN  GIVI 

Department  of  Mechanical  and  Aerospace  Engineering 
University  at  Buffalo,  SUNY 
Buffalo,  NY  14260-4400 


Abstract.  An  overview  is  presented  of  recent  developments  and  contribu- 
tions in  large  eddy  simulation  (LES)  of  turbulent  reactive  flows.  The  foun- 
dation of  some  of  the  recently  proposed  subgrid  scale  (SGS)  closures  for 
such  simulations  is  presented,  along  with  a discussion  of  their  capabilities 
and  limitations.  The  scope  of  the  review  is  limited  to  physical  modeling. 
In  doing  so,  only  issues  pertaining  to  additional  complexities  caused  by 
chemical  reactions  are  discussed.  That  is,  the  challenges  associated  with 
“general”  LES  of  non-reactive  flows  are  not  considered,  even  though  all 
of  these  challenges  are  indeed  present  (and  in  most  cases  are  a lot  more 
complex)  in  reactive  flows.  It  is  recognized  that  numerical  algorithms  and 
computational  procedures  play  a significant  role  in  (any)  LES.  However, 
this  review  does  not  deal  with  these  issues  except  for  cases  wherein  the 
actual  numerical-computational  methodology  is  directly  coupled  with  the 
procedure  by  which  LES  is  conducted.  The  SGS  closure  based  on  the  re- 
cently developed  “filtered  density  function”  (FDF)  method  is  described  in 
a greater  detail.  This  is  due  to  more  familiarity  of  this  reviewer  with  this 
closure;  it  does  not  imply  that  other  closures  are  less  effective. 


1.  Introduction 

In  the  late  1980’s,  I was  preparing  a review  article  on  large  scale  numerical 
simulations  of  turbulent  reactive  flows.  The  intent  was  to  provide  a sur- 
vey of  the  contributions  made  to  both  direct  numerical  simulation  (DNS) 
and  large  eddy  simulation  (LES).  However,  when  that  article  was  finally 
published  (Givi,  1989),  its  content  was  heavily  biased  towards  DNS.  This 
was  not  intentional,  it  just  reflected  the  state  of  progress  on  LES  of  turbu- 
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lent  combustion  at  that  time.  But  with  all  of  the  enthusiasm  for  DNS  in 
the  combustion  community,  the  limitations  of  such  simulations  were  well 
recognized  (even  with  the  most  optimistic  predictions  of  growth  in  super- 
computer technology).  It  was  also  clear  that  the  future  of  large  scale  sim- 
ulations of  practical  turbulent  reacting  flows  would  heavily  depend  on  the 
development  of  LES.  Therefore,  it  was  quite  easy  to  predict  that  LES  would 
receive  significant  attention  in  computational  predictions  of  turbulent  re- 
acting flows  in  the  1990’s  and  into  the  next  (present)  century. 

Now,  at  the  time  of  writing  this  article  (Summer  2001),  while  struggling 
to  meet  the  deadline  for  its  submission(l)  I am  not  surprised  by  the  extent 
of  the  contributions  in  developing  subgrid  scale  (SGS)  models  or  by  the 
magnificent  work  on  LES  of  a variety  of  turbulent  reacting  flow  systems. 
In  fact,  I admit  that  the  rate  of  these  developments  has  been  a lot  faster 
than  my  capability  to  absorb,  or  in  some  cases  even  follow,  the  details  of 
the  proposed  methodologies.  In  addition,  the  page-limit  restrictions  under 
which  this  is  being  prepared,  preclude  describing  the  details  of  the  wide 
variety  of  currently  available  closures;  similarly,  citation  of  the  relevant 
references  cannot  even  be  done  exhaustively.  Fortunately,  many  aspects  of 
SGS  closures  and  LES  of  reacting  turbulence  have  recently  been  discussed  in 
several  excellent  tutorial  and  review  articles  (Cook  and  Riley,  1998a;  Candel 
et  ai,  1999;  Bilger,  2000;  Branley  and  Jones,  2000;  Menon,  2000;  Peters, 
2000;  Pope,  2000;  Luo,  2001;  Poinsot  and  Veynante,  2001).  Therefore,  in 
the  present  review  I concentrate  on  some  of  the  major  issues  related  to  my 
area  of  research  within  this  field. 

2.  Starting  Equations 

Large  eddy  simulation  involves  the  use  of  the  spatial  filtering  operation 
(Sagaut,  2001) 


/+oo 

Q(x,,t)G(x,,x)dx',  (1) 

-OO 

where  Q denotes  the  filter  function  of  width  Ac,  and  (Q(x,t))e  represents 
the  filtered  value  of  the  transport  variable  Q(x,  t).  In  variable  density  flows 
it  is  convenient  to  consider  the  Favre  filtered  quantity, 

(Q(x,t))L=(pQ)e/(p)e-  (2) 

We  consider  spatially  & temporally  invariant  and  localized  filter  functions, 
Q(x',x)  = G(x' — x)  with  the  properties  G(x)  = G(— x),  and  G(x)dx  = 

1.  Moreover,  we  only  consider  “positive”  filter  functions  for  which  all  the 
moments  xmG{x)dx  exist  for  m>  0. 
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To  set  the  framework,  we  consider  the  transport  equations  of  chemically 
reacting  flows.  To  isolate  the  effects  of  chemical  reaction  in  the  simplest 
way,  we  consider  single-phase  (gaseous)  combustion  in  a low  Mach  number 
flow  with  negligible  radiative  heat  transfer  and  viscous  dissipation.  We  also 
assume  that  Newton’s  law  of  viscosity,  Fourier’s  law  of  heat  conduction  and 
Fick’s  law  of  mass  diffusion  are  applicable.  Therefore,  the  primary  transport 
variables  are  the  density  p,  the  velocity  vector  i = 1,2,3  along  the 
Xi  direction,  the  pressure  p,  the  species’  mass  fractions  Ya,  and  the  total 
specific  enthalpy  h.  All  of  the  mass  fractions  and  the  enthalpy  are  grouped 
into  the  scalar  array  <£(x,t)  = [</q,  fa, . . = [Fi,  Y2, . . . , Yfje,  h]  of  size 

a — Ns  + 1 where  Ns  denotes  the  total  number  of  species.  Application  of 
the  filtering  operation  to  the  equations  of  continuity,  momentum,  enthalpy 
(energy)  and  species  mass  fraction  equations  gives 

^ + ^F  = 0.  (3) 

at  oxi 

d{p)t{uj)L  d{p)t{uj)L{uj)L  d(p)(  d(Tij)t  dTjj 

dt  dxi  dxj  dx,i  dxi  ’ ' ' 

d(p)t{<t>a)L  d(p)((ui)L(<Pa)L  d{Jf)e  dMf 

— Tt — + W, ~ Ihi  &T  + (5) 

where  t represents  time,  and  the  filtered  reaction  source  terms  are  denoted 
by  ( pSa)e  - ( p)e{Sa)i . The  viscous  stress  tensor  and  the  mass/heat  fluxes 
are  denoted  by  r,-j,  and  Jf , respectively.  At  low  Mach  numbers  and  heat 
release  rates,  by  neglecting  the  viscous  dissipation  and  thermal  radiation 
the  source  terms  in  the  enthalpy  equation  can  be  assumed  to  be  negligible. 
Thus,  Sa  = Sa(<t>).  The  terms  TtJ  = (p)e((uluJ)L  - (ui)L(uj)L)  and  Mf  = 
(p) (((ui4>a) l — {ui)i{(f)a)i)  denote  the  SGS  stress  and  the  SGS  mass  flux, 
respectively.  Equations  (3)-(5)  are  coupled  through  the  equation  of  state. 

3.  Closure  Methodologies 

For  non-reacting  flows  the  SGS  closure  is  associated  with  Tij  and 
(Canuto,  1994;  Ciofalo,  1994;  Lesieur  and  Metals,  1996).  In  reacting  flows, 
an  additional  model  is  required  for  the  filtered  reaction  rate  ( S0)l ■ This 
modeling  is  the  subject  of  primary  concern  in  this  review. 

One  of  the  first  contributions  in  LES  of  reactive  flows,  similar  to  that  in 
LES  of  non-reacting  flows,  was  made  in  atmospheric  sciences  (Schumann, 
1989).  In  this  work,  the  effects  of  SGS  scalar  fluctuations  (as  appear  in  the 
chemical  source  term)  are  assumed  negligible,  i.e.  (Sa(<f>))L  ~ .S'«((</>)l). 
This  assumption  is  compatible  with  that  made  in  some  of  the  more  recent 
contributions  (Boris  et  al. , 1992;  Fureby  and  Grinstein,  1999),  in  which  it 
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is  argued  that  all  of  the  essential  SGS  contributions  are  included  in  the 
numerical  discretization  procedure. 

Modeling  of  the  scalar  fluctuations  has  been  the  subject  of  broad  in- 
vestigations in  Reynolds  averaged  simulations  (RAS)  for  over  five  decades, 
resulting  in  a variety  of  closure  strategies  (Libby  and  Williams,  1980;  Libby 
and  Williams,  1994).  Within  the  past  10  years  or  so,  almost  all  of  these 
closures  have  been  considered  for  LES.  Examples:  the  eddy-break  up  mod- 
els (Fureby  and  Lofstrom,  1994;  Candel  et  al.,  1999),  moment  methods 
(Frankel  et  al.,  1993),  the  flamelet  concept  (Cook  et  al.,  1997;  Cook  and 
Riley,  1998b;  De  Bruyn  Kops  et  al.,  1998;  DesJardin  and  Frankel,  1998; 
DesJardin  and  Frankel,  1999;  Pitsch  and  Steiner,  2000;  Ladeinde  et  al., 
2001),  the  linear  eddy  model  (LEM)  (McMurtry  et  al.,  1992;  Menon  and 
Calhoon,  1996;  Kim  et  al.,  1999;  Menon,  2000),  the  conditional  moment 
method  (CMM)  (Bushe  and  Steiner,  1999;  Steiner  and  Bushe,  2001),  and 
many  others  (Sykes  et  al.,  1992;  Galperin  and  Orszag,  1993;  Smith  and 
Menon,  1996;  Im  et  al.,  1997;  McGrattan  et  al.,  1998;  Thibaut  and  Candel, 
1998;  Battaglia  et  al.,  2000;  Collin  et  al.,  2000).  In  addition,  several  of  the 
closures  previously  developed  for  LES  of  non-reacting  flows,  have  been  ex- 
tended for  use  in  reacting  flow  simulations  (DesJardin  and  Frankel,  1998; 
Jaberi  and  James,  1998). 

The  probability  density  function  (PDF)  methods  have  proven  particu- 
larly useful  in  RAS  (O’Brien,  1980;  Pope,  1985;  Dopazo,  1994;  Fox,  1996; 
Pope,  2000).  The  systematic  approach  for  determining  the  PDF  is  by  means 
of  solving  its  transport  equation.  An  alternative  approach  is  based  on 
assumed  methods  in  which  the  shape  of  the  PDF  is  specified  a priori. 
This  has  been  pursued  in  several  studies  in  most  of  which  it  is  assumed 
that  the  thermo-chemical  variables  depend  only  on  the  mixture  fraction, 
e.g.  infinitely  fast  reaction,  equilibrium  chemistry.  Therefore,  the  PDF  is 
univariate  (Madnia  and  Givi,  1993;  Cook  and  Riley,  1994;  Reveillon  and 
Vervisch,  1996;  Branley  and  Jones,  1997;  Jimenez  et  al.,  1997;  Mathey  and 
Chollet,  1997;  DesJardin  and  Frankel,  1998;  DesJardin  and  Frankel,  1999; 
Forkel  and  Janicka,  2000;  Kempf  et  al.,  2000).  For  LES  of  non-equilibrium 
reactive  flows,  it  is  necessary  to  assume  the  joint  PDF  of  multi-scalars 
(Frankel  et  al.,  1993).  Consistent  with  popular  methods  of  generating  uni- 
variate (Leemis,  1986)  and  multivariate  (Johnson  and  Kotz,  1972)  distribu- 
tions, all  of  the  assumed  SGS  scalar  PDFs  in  the  contributions  cited  above 
are  based  on  the  first  and  the  second  order  moments.  The  PDFs  generated 
in  this  way  offer  sufficient  flexibility  and  are  affordable  for  large  scale  sim- 
ulations. However,  it  is  now  well  understood  that  the  “true”  PDF  strongly 
depends  on  the  actual  physics  of  mixing  in  a given  flow  condition  (Jaberi 
et  al.,  1996).  Therefore,  there  is  a need  to  determine  such  PDFs  in  a more 
systematic  manner. 
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The  “filtered  density  function”  (FDF)  methodology  introduced  by  Pope 
(1990)  provides  the  framework  for  fundamental  developments  of  the  PDF 
based  SGS  closures.  This  method  provides  a means  of  determining  the  PDF 
from  its  own  transport  equation.  For  the  scalars’  array  0(x,  t)  the  FDF, 
denoted  by  Pl,  is  defined  as  (Pope,  1990) 

/ + CO 

c [0,  0(x',  t )]  G(x'  - x)dx',  (6) 

-OO 


a 

C [Vb  <£(x,  *)]  = n (7) 

Q!=l 

where  6 denotes  the  delta  function  and  0 denotes  the  composition  domain 
of  the  scalar  array.  The  term  Q[<f>  - 0(x,  <)]  is  the  “fine-grained”  density 
(Lundgren,  1967;  O’Brien,  1980;  Pope,  1985;  Dopazo,  1994).  In  variable 
density  flows,  it  is  convenient  to  consider  the  “filtered  mass  density  func- 
tion” (FMDF),  denoted  by  F h,  as 

/ + CO 

p(x',  t)C  [0, 0(x',  <)]  G(x!  - x)dx\  (8) 

-OO 


The  integral  property  of  the  FDF  and  FMDF  is  such  that 

/+oo  r+oo 

Pi(0;  X,  f)d0  = 1,  / Fl(0;  x,  t)rf0  = (p{x,  t))t 

-oo  J — OO 


(9) 


For  further  discussions,  it  is  useful  to  define  the  mass  weighted  conditional 
filtered  mean  of  the  variable  Q(x,  t), 


<Q(x,i)  |0)*  = 


/_  ” P(x',  t)Q(x',  t)C  [0, 0(X;,  t)}  G(x'  - x)rfx' 
Fl(0;x,<) 


(10) 


Therefore,  when  Q can  be  completely  described  by  the  compositional  vari- 
able, i.e.  Q(x,t ) = <5(0(x,  /)),  we  have  (Q(x,t) |0)r  = Q(0).  Also, 


/+oo 

(Q(x,t)|0)fFL(0;x,t)d0  = (p(x,t))e{Q(x,t))L.  (11) 

-OO 


The  transport  equation  for  F^(0;  x,  t)  is  obtained  by  multiplying  the  trans- 
port equation  for  the  fine  grained  density  by  the  filter  function  G(x/  — x) 
and  integrating  over  x'  space  (Gao  and  O’Brien,  1993;  Colucci  et  al,  1998; 
Reveillon  and  Vervisch,  1998;  Jaberi  et  oi.,  1999;  Jaberi,  1999;  Zhou  and 
Pereira,  2000;  Tong,  2001), 
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dFL(ij>-,  x,  t)  d[(uj(x,t)\ip)eFL(^-,-x,t)\ 
dt  dx{ 


d[Sa(il>)FL(iJy,x,t)} 

di>a 


The  first  term  on  the  RHS  is  due  to  chemical  reaction  and  is  in  a closed 
form.  This  demonstrates  the  primary  advantage  of  the  FDF  methodology. 
However,  the  SGS  convection  (the  second  term  on  the  LHS)  and  SGS  mix- 
ing (the  second  term  on  the  RHS)  must  be  modeled.  One  of  the  most  chal- 
lenging issues  in  FDF  is  associated  with  closure  of  the  mixing  term.  This 
has  been  the  subject  of  broad  investigations  in  PDF  modeling  (Pope,  1985; 
Pope,  2000).  In  Eq.  (12)  the  effects  of  mixing  are  displayed  through  the 
“conditional  expected  diffusion”  of  the  scalars,  but  can  also  be  represented 
in  the  form  of  the  “conditional  expected  dissipation”  (O’Brien,  1980;  Pope, 
1985).  The  closure  for  this  can  be  via  any  of  the  ones  currently  in  use  in 
PDF  methods  (Pope,  2000).  In  the  absence  of  a clearly  superior  model,  the 
linear  mean  square  estimation  (LMSE)  model  (O’Brien,  1980)  has  been 
used  in  almost  all  of  previous  LES  based  on  FDF  (Colucci  et  al,  1998; 
Jaberi  et  al.,  1999;  Garrick  et  al.,  1999;  James  and  Jaberi,  2000;  Zhou  and 
Pereira,  2000).  With  Jf  = this  model  is 


d_  \ d(FL/(P)ey 

dx{  . dxi 

r\ 

+ ~ (<Ml)Fl]  ,(13) 


where  f2m(x,f)  is  the  “frequency”  of  mixing  within  the  subgrid  and  must 
be  modeled.  The  convective  term  can  be  modeled  as 


Ui\i>)eFL  = ( Ui)LFL  - 

where  7 1 is  the  SGS  diffusion  coefficient  and  must  be  specified.  Equation 
(14)  is  in  accord  with  that  often  used  in  conventional  LES  (Moin  et  al., 
1991;  Canuto,  1994;  Ciofalo,  1994;  Lesieur  and  Metais,  1996).  With  this 
formulation,  obviously  the  resolved  hydrodynamic  field  must  be  determined 
by  other  means.  This  problem  can  be  circumvented  by  considering  the  joint 
velocity-scalar  FMDF, 


/+00 

p(x',  t)Z  [v,  u(x',  t),  V>,  0(x',  f)]  G(x'  - x)dx',  (15) 

-OO 
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£[v,u(x,<),tf,^(x,f)]  = J]  0]  11  ^at(^G  0]?  (Ib) 

fc=l  a=l 

where  v denotes  the  composition  domain  of  the  random  velocity  vector,  and 
£[v,  u(x,  f),  x[>,  </>(x,  f)]  is  the  fine-grained  velocity-scalar  density.  Most  re- 
cent work  in  this  regard  consider  the  transport  of  the  velocity  FDF  (VFDF) 
(Gicquel,  2001)  and  the  joint  velocity-scalar  FDF  (VSFDF)  (Drozda,  2001). 
The  operational  procedure  is  similar  to  that  developed  previously  for  PDF 
methods  (Pope,  1985;  Pope,  1994;  Pope,  2000). 

The  closure  problems  as  noted  above  are  not  particular  to  the  FDF; 
all  of  the  other  schemes  require  similar  modelings.  For  example,  in  the 
limit  of  equilibrium  chemistry  all  of  the  statistics  of  the  reacting  fields 
are  related  to  those  of  the  mixture  fraction.  The  FMDF  of  the  mixture 
fraction  can  be  obtained  from  the  solution  of  Eq.  (12)  with  5 = 0.  So, 
there  is  still  a need  for  modeling  of  the  mixing  term.  Even  in  cases  where 
the  FMDF  is  assumed,  its  distribution  is  parameterized  with  the  low  order 
moments  of  the  mixture  fraction.  As  indicated  above,  the  first  two  moments 
are  typically  used  for  this  parameterization.  Therefore,  there  is  a need  for 
closure  of  the  “total  SGS  dissipation”  as  appears  in  the  second  moment 
(SGS  variance)  equation.  Several  means  of  dealing  with  this  closure  problem 
are  available  (Girimaji  and  Zhou,  1996;  Pierce  and  Moin,  1998;  Cook  and 
Bushe,  1999;  Jimenez  et  al.,  2001;  De  Bruyn  Kops  and  Riley,  2001). 

The  above  problem  is  a bit  more  complex  when  the  SGS  chemical  reac- 
tion is  assumed  to  be  in  the  “flamelet”  regime  (Peters,  2000).  In  this  case, 
even  with  the  one-dimensional  flamelet  model,  the  thermo-chemical  vari- 
ables are  parameterized  by  the  mixture  fraction  and  its  rate  of  dissipation 
(Cook  et  al.,  1997;  Cook  and  Riley,  1998b;  De  Bruyn  Kops  et  al. , 1998; 
DesJardin  and  Frankel,  1998;  Cook  and  Riley,  1998b).  Therefore,  there  is  a 
need  for  a priori  specification  of  the  joint  FDF  of  the  mixture  fraction  and 
its  dissipation.  A review  of  different  methods  of  dealing  with  this  issue  is 
available  (Cook  and  Riley,  1998a).  Equation  (12)  with  5 = 0 indicates  that 
there  is  a dependence  between  the  FDF  of  the  mixture  fractions  and  the 
conditional  expected  diffusion  (and  the  conditional  expected  dissipation). 
This  dependency  is  not  considered  in  most  previous  contributions,  but  is 
the  subject  of  current  investigations  (DesJardin  et  al,  2001). 

Modeling  of  the  conditional  expected  dissipation  is  also  required  in  the 
conditional  moment  method  (Bushe  and  Steiner,  1999;  Steiner  and  Bushe, 
2001).  This  issue  has  been  recognized  at  the  early  stages  of  developments  of 
CMM  in  RAS  (Bilger,  2000).  With  this  model,  the  conditional  filtered  mean 
values  of  the  thermo-chemical  variables  (LHS  of  Eq.  (10))  are  obtained 
by  their  modeled  transport  equation.  This  is  obviously  computationally 
less  demanding  that  solving  the  FDF  transport  equation.  But  in  order  to 
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determine  the  actual  filtered  quantities,  the  distribution  of  the  mixture 
fraction  FDF  must  be  specified. 

An  important  issue  in  regard  to  FDF  is  associated  with  the  numeri- 
cal solution  of  its  transport  equation.  The  Lagrangian  Monte  Carlo  scheme 
(Pope,  1985)  has  proven  particularly  useful  for  this  purpose.  In  this  scheme, 
the  FDF  is  represented  via  an  ensemble  of  computational  elements  or  par- 
ticles. Transport  of  these  particles  and  the  change  in  their  properties  are 
modeled  by  a set  of  stochastic  differential  equations  (SDEs)  (Soong,  1973). 
The  diffusion  process  (Gardiner,  1990)  has  proven  effective  for  this  pur- 
pose. The  coefficients  in  the  Langevin  equation  governing  this  process  are 
set  in  such  a way  that  the  resulting  Fokker- Planck  equation  (Risken,  1989) 
is  equivalent  to  the  FDF  transport  equation.  Therefore,  the  Monte  Carlo 
solution  of  the  SDEs  represent  the  solution  of  the  FDF  in  the  probabilistic 
sense.  This  procedure  has  proven  successful  for  simulating  PDF  in  a variety 
of  systems  (Grigoriu,  1995).  However,  one  must  be  careful  in  performing 
stochastic  simulations  in  conjunction  with  modern  CFD  solvers.  Many  of 
the  advanced  discretization  routines  developed  for  solving  deterministic  dif- 
ferential equations  may  not  be  applicable,  or  may  have  to  be  significantly 
modified  to  be  suitable  for  solving  SDEs  (Kloeden  and  Platen,  1995). 

Implementation  of  LEM  is  also  based  on  stochastic  representation  of 
the  flow.  In  its  original  development  in  RAS  (Kerstein,  1988),  the  processes 
of  molecular  diffusion,  chemical  reaction  and  turbulent  convection  are  con- 
sidered separately.  This  is  achieved  by  a reduced  one-dimensional  (linear) 
description  of  the  scalar  field,  which  makes  it  possible  to  resolve  the  flow 
scales  even  for  flows  with  relatively  high  Reynolds,  Schmidt  and  Damkohler 
numbers.  The  interpretation  of  the  one-dimensional  domain  is  dependent 
on  the  particular  flow  under  consideration.  In  this  way,  the  processes  of 
molecular  diffusion  and  chemical  reaction  are  taken  into  account  exactly, 
but  the  effects  of  convection  are  modeled.  This  is  achieved  by  “random 
rearrangement”  (or  stirring)  events  in  such  a way  that  the  displacements 
of  fluid  elements  result  in  a diffusivity  equal  to  the  “turbulent  diffusivity.” 
For  LES,  this  procedure  is  followed  within  each  of  the  computational  cells, 
and  stirring  is  performed  to  yield  the  desired  SGS  diffusivity.  Menon  and 
colleagues  have  made  extensive  use  of  LEM  for  LES  of  a wide  variety  of 
reacting  flows.  A recent  review  is  available  (Menon,  2000). 
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